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Abstract 


The quantum Liouville equation in the Wigner representation is 
solved numerically by using Monte Carlo methods. For incremental 
time steps, the propagation is implemented as a classical evolution in 
phase space modified by a quantum correction. The correction, which 
is a momentum jump function, is simulated in the quasi- classical 
approximation via a stochastic process. In this paper the technique, 
which is developed and validated in two- and three-dimensional mo- 
mentum space, extends an earlier one- dimensional work. Also, by 
developing a new algorithm, the application to bound state motion in 
an anharmonic quartic potential shows better agreement with exact 
solutions in two-dimensional phase space. 


I. Introduction 

Nuclear interaction theory is formulated in the language of quantum mechanics, and hence 
the development of general methods of solutions to quantum dynamics will provide useful tools 
for application to a large class of problems in nuclear many-body theory. Different approaches 
exist to the formulation of this complex problem, and attempts toward solutions at various 
levels of approximations are ongoing. The time-independent approach based on the Lippmann- 
Schwinger equation, for instance, is useful for describing systems with well-defined incident initial 
states (ref. 1). Similarly, the time-independent classical transport theory provides a method for 
calculating the fluence of particles as a superposition of sharply defined incident states under 
steady-state conditions (ref. 2). In real-life situations in deep-space radiations, for example, 
sporadic bursts of radiation may be encountered during which interactions and scattering proceed 
as fast transient events. This is in contrast with the slow variation of background space radiation 
for which a time-independent approach to study the effects of radiation on spacecraft is indicated. 
With this report we begin the development of a practical numerical code that is designed for the 
study of time-dependent nuclear scattering and interaction for transient thermodynamic wide 
spectrum radiation that is aimed toward application to the NASA radiation protection program 
for space travelers. 

The density operator formalism of quantum dynamics (ref. 3) provides a suitable framework 
for the study of thermodynamic systems. In the Wigner representation (refs. 4-7), the dynamic 
equation of the density operator, given by the quantum Liouville equation, is transformed into 
ordinary functions and operators in phase-space coordinates. In a series expansion in powers of 
Planck’s constant h, the equation then provides an intuitively appealing reduction to the classical 
Liouville equation in the classical limit. Also, the more familiar equations appearing in nuclear 
scattering and heavy-ion collision theory, such as the hydrodynamic equations (see refs. 5 and 8)- 
and the Boltzmann- Vlasov equations, may be extracted from the Wigner formalism. Because 
many of the cross sections used in the space program are derived from Monte Carlo simulation 
of the classical Boltzmann transport theory (ref. 9), the quantum correction to classical theory 
is of interest to NASA. 

In this paper Monte Carlo methods are applied to solve the quantum Liouville equation in 
the Wigner representation (refs. 10 and 11). The equations are in a noncovariant form and apply 
to single-particle dynamics only. The time evolution is treated as a stochastic process, as seen 
in references 7 and 10-12. In an effort to simplify the problem, only first-order quantum effects 
are considered; and in this approximation the solution is applicable to quasi-classical systems 
(refs. 11, 13, and 14) that exhibit smoothly varying momentum distribution typical of highly 
mixed thermodynamic systems. In general, however, the first-order quantum correction may 


not be sufficient and may, in some instances, even require the entire series summation (refs. 15 
and 16). For the scattering of a highly collimated beam, for example, higher order terms become 
increasingly significant. Therefore, the method pursued in this work will hopefully complement 
the other approaches mentioned earlier. 

A generalized Monte Carlo method was introduced in references 10 and 11. This paper 
extends that work to two and three dimensions, and a new algorithm is developed that gives 
improved results for the application considered. 

In section II, quantum dynamics in the Wigner representation is reviewed and the stochastic 
techniques are developed. In section III, the technique is validated independently of the 
classical motion by comparing it with analytic solutions in the one-, two-, and three-dimensional 
momentum space. In section IV, an application to bound state motion within an anharmonic 
quart ic potential in two-dimensional phase space is considered and the algorithm is discussed. 
In section V, the results and a discussion are presented, and in section VI, some concluding 
remarks and future applications are briefly indicated. 

II. Theory 

Quantum Liouville Equation (QLE) 

The density operator p of a quantum thermodynamical system is given by 

P ~ ^ ^ Pinl^m X | (1) 

m 

where P m is the probability for an ensemble element to be in eigenstate | >. The time 
evolution of p is the quantum Liouville equation, 



where t denotes time, H is the Hamiltonian, and h = 1. Equation (2) has the formal solution 

p(t)=e~ mt p( 0)e^ (3) 

Because the components of H are usually noncommutative, this form is difficult to solve in 
practice. An intuitively appealing solution can be obtained by taking the Wigner transform of 
the QLE, which provides a series expansion in ft and reduces to the classical Liouville equation 
in the classical limit, h 0. 


Wigner Representation of QLE 

A few basic properties of the Wigner transform are now reviewed. The Wigner transform of 
an operator, O is defined by 


Or 


(x,p,f) = J 


dye 


*P T < x - -y 

T 


0(t) 


X+ 2 y> 


( 4 ) 


which is a simultaneous representation in both position coordinates x and momentum coordi- 
nates p. The Wigner transform for the density operator p is 


f°° . 1 1 

fw (x, P, t) = / dy e ip y < x - -y |p (<)| x + -y > 
7—oo ^ ^ 


( 5 ) 
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and is called the Wigner distribution function (generally defined with a normalization fac- 
tor (2tt) -3 ). As an example, the Wigner transform of the density operator corresponding to 
a minimum wave packet defined by 


^(x) 


(27TCT 2 ) 


2 ^ 3/2 


exp 


(x - x 0 )' 

-rpo-x — 2 r - 


( 6 ) 


is given by 


fw (x, p, t) = J dy e ,p ' y xp* (x - ip (x + lyj 


= 2 exp 


-2cr 2 (p - po) 2 


(x - x 0 ) 2 
2 ct 2 


(7) 


The Wigner function has many analogs with the classical distribution function. 

For example, 

( 2 tt )' 3 J dp f w (x, p, t) = < X |p| x > 

(8) 

( 2 tt ) -3 J dx f w (x, p, t) = < p |p| p > 

(9) 

( 27 r)“ 3 J dx dp f w (x, p, t) = 1 

(10) 

and the expectation value of an observable O is given by 


< 0 (£) > = ( 2 tt )“ 3 [ dx dp O w (x, p, t) f w (x, p, t) 

(11) 


However, even though / u ,(x. p, t) is real (that is, f* = f w ) it cannot strictly be a distribution 
function because it can have negative values, and therefore the Wigner function should at most 
be considered as an auxiliary function that is useful for calculating thermodynamic averages. 

The Wigner transform of the quantum Liouville equation becomes 

^(x, p ,t) = —2H W sin / w (x,p,t) (12) 

in which H w is the Wigner transformed Hamiltonian and A is the Poisson bracket operator given 
as _ v 

A = V p • V x - V x • V p (13) 

where the arrows indicate the direction of action of the operator. Expanding the sine term gives 
the series expansion 

(x,p, t) = ^-H W A + — H W A 3 - fw (x, p, t) - (C c + C q ) f w (x, p, t) (14) 

where C c — H W A is the classical Liouville operator and - C q (which is equal to all higher order 
terms) is the quantum operator. The solution to equation (14) is given by 

fw (x, p, t) = e-( £c+ ^)7ic (x, P, 0) (15) 
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For small increments of time, equation (15) becomes 


fw (x, p, At) = e £ " A 'e £cA % (x, p, 0) + O (At 2 ) (16) 

Hence, infinitesimal time motion can be described in terms of successive classical and quantum 
evolutions in which the classical operator transforms the function to 

five (x, P, At) = e~ CcAt f w (x, p, 0) (17) 

and the quantum operator acts on f wc , thus giving 

fw (x, p, At) = e~ L< i At f wc (x, p. At) (18) 

These expressions are difficult to evaluate analytically for arbitrary functions. Hence, Monte 
Carlo methods are applied with the advantage that the only analytic evaluation required is that 
for the action of the operators on a delta function. Explicit expressions for the operators C c 
and Lq for a Hamiltonian operator of the form 


with Wigner transform 

Hw (x, p) = J^ + F(x) 

(where x and p are now variables and not operators) are obtained as 


(19) 


( 20 ) 


Cc = 





( 21 ) 


C n 


24 


H(x) 



( 22 ) 


The action of C c on the delta function is standard. The action of Cq on the delta function is 
now evaluated explicitly in the quasi-classical limit for a central potential V(r) with r = |x - x,;| 
and, 



(23) 


Expanding Cq in terms of the radial component pg and the perpendicular components pi 
and p 2 of the momentum (appendix A) gives 


where 


with 


£q — £q\ + £<72 

(24) 

£</i = &pq + &T dpo 

(25) 

c n = y + a T 0 P{) dj 2 

(26) 

1 d 3 V 


° L ” 24 c?r 3 

(27) 
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Consider first the action of C qi on the delta function. (Similar arguments hold for C qr ) It 
acts on the pq and p\ components only, thus giving 

J (P - Pi) = e~ Cqx *6 (P - Pi) 

= 6 (P2 ~ P2i) e ~ Cqx ^ (PO ~ POi) t (Pi ~ Pit) ( 29 ) 


A change of variables to vq = pq + 7Pi and v\ = po — *yp\ reduces the expression to a product 
of one-dimensional forms. The operator £ qi transforms into 



C qx — a dy Q + a 

(30) 

where 

a = ^ + a T 7 2 = 2 a L 

(31) 

with 

(Za L \^ 

7 V 2 a T ) 

(32) 

and equation (29) transforms into 



J (P - Pi) -»■ 276 (P2 - P2i) e atdv ° <5 («0 - VQi) e atd[, i d («i - Vu) 

(33) 


where is the Jacobian of the transformation for the delta functions. 

The expressions to be evaluated are of the typical one-dimensional form with at replaced 
by a. Thus, 

Ai (a; p - Pi) = e~ ad 'v 6 (p - p t ) 

1 r co o . 

= — / dye iay + ipy (34) 

J —(X) 

where Ai is recognized as the Airy function. Depending on the sign of a, the function decreases 
exponentially along one direction and is oscillatory along the other with a slow decay in amplitude 
and increasing frequency. 

Monte Carlo Method 

In a Monte Carlo procedure a sample set of test points is selected to represent the initial 
positive valued function. Thus, 


3 tV 

fw (x, p, 0) « °i 6 ( x - x *) 6 (p - Pi) ( 35 ) 

i = 1 

where = 1 is the sign of the test point because, as noted above, /^(x, p,£) may be positive 
or negative. The classical propagation is a canonical contact transformation that transports the 
delta functions to new positions along deterministic trajectories so that 

fwc (x, p, At) = &i6 (x - x d ) S (p - Pd) (36) 

i 
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where x c ^ and p (1 are evaluated via Hamilton’s equations of motion 


dxj OH 
clt dp i 

dpj _ dH 
dt dxj 

To implement the quantum correction as given by equation (34), a strong damping for the Airy 
functions is useful. (Details for only the one-dimensional quantum jump function are discussed 
here. For higher dimensions see section III.) Now the phase space assumes a graininess due to 
the delta function representation; that is, the larger the number of representative points, the 
finer the grain structure. For a coarse-grained analysis of the quantum correction, note that the 
increasing rapidity of the oscillation of the Airy function at large momentum distances implies 
a net cancellation. Hence, to speed simulation, a grain size is introduced into the 6(p — p t ) term 
to produce a faster damping rate for the function. This is achieved by approximating the delta 
function by a narrow- width Gaussian function, which modifies the Airy function to 

Ja (a;p- Pi) = e~ atd P8 a ( p - Pi ) (38) 


( 37 ) 


where 8 a (p — pi) is the Gaussian function of width a. The expressions for the modified Airy 
functions J a are given in appendix B. 

The corresponding quantum jump function is defined as 


J a (a;p - Pi) = Ja (a; P - Pi) - 8a {p - Pi) 


(39) 


Figure 1(a) illustrates a typical Gaussian-modified Airy function, and figure 1(b) illustrates the 
corresponding quantum jump function J tt (a;p )> which is shown as “«/” in the figures. 

The jump function is implemented via a stochastic simulation. To this end, let J± correspond 
to the positive and negative segments of the function J a . Partial integration easily shows that, 

J J a (a; p - p^ dp = 0 (40) 

which indicates that the areas under the positive and negative segments are equal. Defining the 
area A gives 

A = J \J±\ dp (41) 

and rewriting equation (39) by using equations (40) and (41) gives 


Ja { a iP ~ Pi) = A 


\J± \ 

A 


\J-\- 

A ; 


= A[F + -F-} 


(42) 


which defines the jump “probability” functions as F±(a;p - p^ = | J-t| /A. 

The stochastic method is based on the following probabilistic interpretation. For the two 
random variables X and T, the joint probability P{X, Y) is given by 


P{X,Y) = P(X\Y)P{Y) 


(43) 
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where P(Y) is the probability for the event Y and P{X\Y) is the conditional probability for 
the event X, provided that event Y has occurred. Compare equation (43) with equation (42). 
If A < 1, interpret P{Y) - A as the probability for the quantum event, or as the creation 
probability. In other words, only the test points selected randomly with probability P(Y) 
undergo quantum events during each time interval. (That is, if A > 1, let A = n + A, where n is 
an integer, and A < 1. Then, the test point will undergo quantum events n times and P(Y ) A 
will determine whether an additional quantum event should take place.) Generally, a is small 
enough (see fig. 2) to ensure that A < 1 so that, at most, one quantum event occurs per time 
step. 

The conditional probability P(X |T) = F+ — F— represents the momentum jump probability 
corresponding to the random variable X = p. A pair of values A p± is selected randomly by 
using the cumulative distributions for F±. In the Monte Carlo representation, this becomes a 
test pair with coordinates, 

6 (p - ( Pd + A p±)) 8{x- Xd) <7i<7± 

where a± = ± 1 for the positive and negative points. The newly created points are appended to 
the initial set to undergo subsequent classical and quantum motions. If A <SC 1, a factor M is 
introduced to enhance the creation probability to MA, with a normalization factor 1/M for the 
new pairs. 

Clearly, in the absence of the classical motion, the stochastic process is a Markoff process. 
That is, with t n -\ < t n , 

F{p{t n )<Pn\p(t),t<t n -l} = F{p(tn)<Pn\p{tn-l)} 

The jump probability for each test point is thus independent of its past history, and depends 
only on its present location in momentum space, 

III. Validation of Stochastic Quantum Motion 

The validity of the technique developed in the previous section is established by comparing 
stochastic quantum time development in momentum space with analytic solutions. This is easily 
done when the initial function is a Gaussian. 

One-Dimensional Quantum Motion 

The quantum time development for the interval t in the quasi-classical approximation is given 
by 

f(p,t) = e~ atd Pf(p, 0) (44) 

With the initial function given by 

/(p,0) = ^e-» 2 / 2 « 2 (45) 

v27ra 

the analytic solution is the Gaussian modified Airy function given in appendix B. 

For the stochastic evolution, a representative set of points for the initial function is chosen 
as follows: A pair of values (p t , fo ) is selected randomly within a specified boundary for / and p 
such that / lies well within the defined area. The function / varies from 0 to /max = 1/ (\/27ra). 
If /(Pi) < /j, then pi is selected; otherwise, it is discarded. Hence, 

/(p.o) « ^rZX'(p-Pt) ( 46 ) 
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where the test Gaussian functions have width a' with a 1 -C a. By dividing the total time t 
into K discrete time intervals (At = t/K), the time development is written as 

f(p,t) = (e- aM $) K f(p, 0) (47) 

During each time step, statistical test points are selected with probability A (see the discussion 
following eq. (42)) and the new test pairs are created at p l + A p±, where values of Ap± are 
Selected with conditional probability F± which get appended to the main list. The updated list 
is propagated in the subsequent time interval. 

In the actual algorithm, the momentum space is divided into grids and the test points are 
assigned on it. With at — 0.1 and a At = 0.001, 100 time steps are executed. The creation 
probability is enhanced by an arbitrary factor M that is set at M = 10 000/7V. Thus, for 
100 initial test points the creation probability is increased 100 times; that is, the smaller the 
number of initial points, the larger the number of pair creations. Each representative pair for 
J a (aAt;p — Pi) is therefore given by 

jjlHp- ( Pi + Ap + )) -6(p- ( Pi + Ap_))] 

For a density k on the grids, the process is repeated k times. 

Figure 3 compares the results with the analytic solutions for various grid sizes, the initial 
number of test points iV, and for various Gaussian widths a f for the test points. The results 
show good agreement with the analytic solutions and appear to be independent of the variables. 


Two-Dimensional Quantum Motion 

The two-dimensional quantum motion is given by 

f(POiPbt) = e~ c ^f(pQ,pi,oy (48) 

where pq and p\ are the radial and perpendicular components, respectively, and C q is given by 
equation (25) with a^j 2 replaced by a £. The initial function is chosen to_be 


/(P0>P1*0) 


w exp 


zM±dl 

2a 2 


(49) 


To obtain the analytic solution, change the variables to vq = pq + pi, and tq = p 0 — p\. Thus, 





(50) 


which is recognized as a product of one-dimensional forms. The inverse transformation is then 
computed to get the analytic solution. (See fig. 4.) 

For the stochastic evolution, consider the action of C q on a test point during the subinterval At 
given by 

J a (a At: p - p,) = e~ C( i At S (po - PQi ) 6 (pi - p u ) (51) 
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Transforming to variables vq and v\ as before gives 

J a {aAt: p-pi) -> 2 

= 2 \j a ) («0 - «0t) + S a' ( v 0 - v 0i)] Wa' ( U 1 “ U 1 *) + 6 a' ( U 1 ~ V h)] ( 52 ) 


e a Atd ’ v ^S a i (vq - vq, ) 


e aAtd ^S a ,(vi-v u ) 


which to 0(At) gives 

J a (a A t; P - Pi) = 2 [J a , (v 0 - i>oi) ^ (”1 “ ”h) + J a' ( v i “ ^ (**> - t>0i) + fa) - «0i) ^ fa ” v i01 

(53) 

because J Q / is of O(At). The pair selection for each J a t is done as before and the representative 
test pairs are 

6 {v\ — v\i) 6 (uo — (^Oi + A y Q±)) o±Gi + <5 (r>o u 0z) ^ fal ~ fal* ^ v l±)) cr ± (J i 

Note that two pairs are created for each event expressed by the two summations. Transforming 
back to the original coordinates gives the representative test pairs 

6 (po ~ («h + ^)) 6 (pi - (ph + ^)) +*(*>- (m + ^r)) 6 (pi - (pu - ^)) 

Figures 5(a) and 5(b), which show the results for stochastic simulation, compare well with 
figures 4(a) and 4(b), respectively. As before, 100 time steps were executed, and the pair creation 
probability was enhanced by a factor of 20 by using an initial number of 10 000 test points. The 
effect of increasing the width parameter a ' on the simulation is seen by comparing figures 5(a) 
and 5(b) with figures 5(c) and 5(d). The effect of increasing grid size is seen in comparing 
figures 5(a) and 5(b) with figures 5(e) and 5(f). Although a 25-percent increase in width a' has 
little effect on the solution, the use of a 33-percent larger grid lowers the distribution peaks, as 
can be seen when comparing figure 5(f) with figure 5(b). 


Three-Dimensional Quantum Motion 

The three-dimensional quantum motion is given by 

/(p,t)=e-^/(P>0) (54) 


The initial Gaussian function may be written in terms of parallel and perpendicular components. 


Thus, 


/( P»0) 



(55) 


Similarly, from appendix A, 

C q = a i dp Q + ap d Po d p± (56) 

Hence, the analytic solution is similar to the two-dimensional case on a plane defined by pa 
and p’x- For the stochastic time development, consider C q acting on a test point during time 
interval At. Thus, 


J(p- Pi) = e~ £<? 2 At e £ At 6 (po - Poi) & (Pi ~ Ph) s (P2 - P2i) (57) 

where £ qi and C q2 are given by equations (25) and (26), respectively. The sample set generated 
by £ qi and C q2 acting successively on the test point creates four new pairs to O(At). 
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The operator C qi generates two sets of pairs as in the two-dimensional case that can be 
written succinctly as 


2 

6 (P0~ POj) 6 {pi - Plj ) 6 (P2 - P2i) <Ti<Tj± 

j &= 1 

Similarly, C q2 acting on <5(p - pj) generates the set, 

2 

z2 6 (po - POj) <5 (P2 - P2j) <5 (pi - Pli) 0i<Tj± 

j &= 1 

Figure 6 shows the results for the (po, pi) plane. The comparison with analytic solutions (fig. 4) 
is remarkably good even with 10 000 initial test points. Also, by choosing a At = 0.01, only 
10 time steps are required. 

IV. Application in Two-Dimensional Phase Space 

The full quantum motion, namely the classical evolution followed by the quantum jumps, is 
applied to an arbitrary initial state in an anharmonic quartic potential: 

V (x) = i (x 2 + fcx 4 ) (58) 

Note that this potential provides an exact description of the quantum effects within the quasi- 
classical approximation as all higher order terms vanish. The problem is first studied in two- 
dimensional phase space to validate the technique with exact solutions calculable by standard 
numerical techniques. The power of the technique developed herein lies in its direct applicability 
to higher dimensions and to many-body systems. 

The initial Wigner functions are chosen from a class of functions represented by 

/ (x, p, 0) = 2 exp {-/? [(x - x 0 ) 2 -{p- po) 2 ] } (59) 

such that / w = Pf, where the parameter 3 defines arbitrary admixtures of states. The examples 
considered have x 0 = 0, and p 0 = 1. (See fig. 7 for /? = 0.25.) With 0 = 1, the Wigner function 
corresponds to a minimum wavepacket that is a pure state. (See eq. (7).) For 0 < 1 the function 
therefore describes a mixture of states. Obviously, 0 > 1 is not allowed because of the uncertainty 
relations Ax A p < | . 

The algorithm is based on the following complex of procedures using C-language. The initial 
set of test points is assigned to a fine mesh of phase-space grids. A list of structures is constructed, 
each structure containing the data corresponding to the coordinates of the grid, the density, and 
the sign of the test points. Only the nonempty grids form the list. For the classical motion 
with mass m = 1, the coordinate data are updated by using a two-step second-order Runge- 
Kutta method. This computation can be as accurate as desired and does not involve a grid 
approximation. A high degree of accuracy is essential for the classical motion. 

To implement the quantum event, all test points within a particular region in position space 
having all possible momentum values are identified by sorting. (To facilitate sorting, the list 
is constructed at two levels. The first level, which consists of structures for a coarse x-grid, 
forms the main trunk. From each unit on the trunk, a branch containing all the structures that 
fall within that unit are attached. The second-level structures contain the actual data.) The 
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selected set is then allowed to undergo one-dimensional quantum jumps. (See section III.) The 
cumulative distribution for F±(a:p) is tabulated for various values of a. The required value of a 
is computed at the coarse x-grid location via a = [V"'(x) At]/24. The net sum of newly formed 
test points is attached to the main list. The entire x-space is spanned in this manner. 

With low creation probabilities and the annihilation of pairs of opposite signs within the 
assigned grid spacing for quantum motion, the main list does not increase exponentially and 
remains tractable. The initial number of test points ( N ) was taken to be 20000, which formed 
an initial list size of approximately 4000 and grew to a size of approximately 15000 at the end 
of t — 4n. The enhancement factor M was chosen as M = 5 with the grid size (annihilation 
distance) set at approximately 0.3. The test points were given Gaussian width a — 0.4. The 
function is reconstructed at the required time intervals from the test points by using a suitable 
set of orthonormal harmonic oscillator test functions. On a micro VAX-4000 series computer 
(manufactured by Digital Equipment Corporation), the run time for the 0 -tt time segment was 
typically 5 minutes, but for the 0-47T time segment it was approximately 40 minutes because of 
the increasing list size. 

An earlier version of the algorithm was written in PL/I language (ref. 10). One complicated 
feature of the algorithm was the task of keeping track of the four nearest neighbors of a moving 
sample test point in order to facilitate sorting and annihilations of the newly created pairs with 
their nearest neighbors having opposite signs. The algorithm developed here has proven to be 
faster and more accurate. 

V. Results and Discussion 

Snapshots of the motion at time intervals in units of 7r are shown in figures 8-14 for various 
initial Wigner functions and for various strengths of the potential. Each time unit is subdivided 
into 30 time steps. The results are compared both with the exact solution calculated by standard 
numerical techniques (ref. 10) and with the solutions of the classical Liouville equation. 

The following observations can be made regarding the classical motion versus the quantum 
motion. For the classical motion, the volume of phase space occupied by the system (an integral 
invariant of Poincare) remains constant (ref. 17), but it streams out into all phase-space regions 
allowed by energy conservation, with the occupied phase-space region developing whorls and 
tendrils. (See fig. 14.) After long intervals of time, this spread gives the appearance of a uniform 
distribution over a coarse grid, although finer grids would reveal the fine detail of the contour 
levels as they are the classical solutions shown in figures 8-14 (part (a)). For the quantum 
motion, however, the system maintains a cohesiveness as the unit oscillates within the potential 
well. This cohesion is the result of quantum interference effects arising from the oscillations of 
the Airy functions, thus causing cancellations and reinforcements over the classical motion. 

Quantitative differences for the pure state (0 = 1) and the mixed state (3 < 1) quantum 
motions are also evident. For the pure state motion, the maximum height of the Wigner function 
is observed to remain unchanged. However, the mixed-state motion shows a “quantum focusing” 
effect as the Wigner function peaks beyond its initial maximum. Clearly, classical motion does 
not allow for such effects resulting from the Liouville theorem, which states that the density 
of systems in the neighborhood of some given system in phase space remains constant in time 
(ref. 17). 

Finally, as an example of computation of an observable quantity, the averages of x and p are 
shown in figure 15, where 

< x > = (27 t) _1 [ dxdpfw (x,p, t) x (60) 
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<P>=(2tt) 1 j dxdp f w {x,p,t)p (61) 

The averages are plotted both for the purely classical and the full quantum motion. For the 
classical motion, the system distributes uniformly around the equilibrium point, consistent with 
energy conservation, and the first moments of the distribution approach zero at late times. For 
the quantum motion, however, these moments are oscillatory with finite amplitude, an indication 
of a preservation of structural unity over long intervals of time. 

Statistical fluctuations are inherent in any Monte Carlo simulation. By increasing the number 
of initial test points, these fluctuations can be made negligible and a single computer run then 
becomes sufficient for accuracy. In conclusion, the method pursued in this work shows great 
promise for application to multidimensional problems in which other numerical procedures may 
prove to be difficult. 

VI. Concluding Remarks 

The quantum Liouville equation in the Wigner representation is solved numerically by 
using Monte Carlo methods. For incremental time steps, the propagation is implemented as 
a classical evolution in phase space modified by a quantum correction. The correction, which is 
a momentum jump function, is simulated in the quasi-classical approximation via a stochastic 
process. In this paper the technique, which is developed and validated in two- and three- 
dimensional momentum space, extends an earlier one-dimensional work. Also, by developing a 
new algorithm, the application to bound state motion in an anharmonic quartic potential shows 
better agreement with exact solutions in two-dimensional phase space. 

Work is well under way toward the development of a code to a six-dimensional case for 
application to potential scattering problems and low-energy barrier penetration. Future work 
will involve extensions to few-body scattering and the inclusion of quantum statistics to account 
for the Pauli blocking effects of spin one-half fermion systems. These are long-term projects, 
but a beginning has been made. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
November 23, 1993 
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Appendix A 


Quantum Evaluation Operator 

In this appendix we expand C q in terms of the radial component po and the perpendicular 
components p\ and p 2 of the momentum. For this we evaluate 

Q = V(r) ( V x ■ V p ) 8 (p) = V (r) J ( V x ■ y) e*™dy 

in terms of parallel (0) and perpendicular (_L) components to get 

V x = ^e r 9 r , e x 

y = (e r yo,e x y_L) 


Using the relations 


d x yo = y± d±y± = -3/0 


gives 


(Vx • y) V (r) = V'yo 


V 7 


(V x -y) z F(r) = E"yo + yyi 

(V x • y) 3 V (r) = V"'yl + 3d r (y) y 0 y\ 


Q (2? r) 3 / 


v'"yl + 3 dr ( — J yoy x 


V 1 


JP 


y dy = 


V'”dl Q +3d r 


(?) 


Opo ^Px 


<5(p) 


Hence, 


Cq ~ 24 


F"'a 3 0 + 3 dr 


(?) 


dpo i 
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Appendix B 

Expressions for Modified Airy Functions J a 

Expressions for the damped Airy functions are obtained in this appendix. The expression to 
be evaluated is 

Ja (a; P ) = e~ ad P 6 a (p) 

where 

v27ra 

The results are presented here. For details, see reference 10 where the evaluation is done by 
using the method of steepest descent (ref. 18). 


Series Expansion for a 1 


~-ad% 


Using the series expansion for e p and applying the Rodrigues formula gives 


where H is the Hermite polynomial. 

Asymptotic Evaluation for p — ► oo 

If the aforementioned expression is rewritten using the integral representation for a Gaussian 
function 


6 a (p) = 


1 


V2 


i: 


7TQ J — oo 


dy exp 


u/2py\ 
a - y 2 J 


the result is 


Ja (a;p) = Re 


^ fVxp 
Jo 


KCt Jo 


v y3 ~ W\ ±iv , 


where 


, V2\ , (y/T 

p -brJ p a = hr 


The integral is evaluated by the method of steepest descent in the complex plane. If we define 

Jl 


f (z) = ia'z 3 


VWl 


± iz 


where 2 is a complex number, the saddle points occur at 


z 0 


j= 1 ± (l + 3a p) 1 ^ 2 ] 

V V] J 


The integral is evaluated independently along different paths for two cases, and the resulting 
expressions are given as follows: 
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For 1 + 3 a'p' > 0, 


Ja (a; p) = Re— 3 = — exp 
lira 


J |3/2 


/ (^o) ^2 


00 (ia'^rpn+l)/^ 

^n!|l + 3aVp +1 )/ 4 


For 1 + 3 a'p' < 0, 


\/2 

Ja (a; p) = Re — exp 

'KOl 


|p'I 3/2 /w + t E 


(ia t e s "'i , \ n T [(3n + 1) /2] 


4 1 newn=0 »'■ U + 3«VI <3 ” +1)/4 


In the region (1 + 3a V) « 0 with 3a'p r < 0, the resulting expressions are given as follows: 
For 1 + 3a 1 p' > 0, 


s/2 

Ja ( a;p ) = Re- — exp 

oTTOt 


»f / 2 fw+% E 


[- (1 +3aV) 1/2 e" /3 ]"r[(2n+ l)/3] 


n~ 0 


n\a 


, ,(2«+l)/3 


For 1 + 3aV < 0, 


%/2 

Ja (a;p) = Re— exp 


x < exp 


oo 

|p'| 3/2 /(2o)] E 

71=0 

* ( n + - 1 - I + exp 


(1 + 3 a V) n/2 r [(2 n + 1) /3] 


n\a 


| /(2n+l)/3 


0 f]+“p[ i (»+5)’ r ]} 
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(a) N = 100; a' = 0.3. 



P 


(b) N = 10 000; a' = 0.3. 

Figure 3. Stochastic evolution of jump function J«(a;p) after 100 time steps compared with 
analytic solution (solid line) at a = 0.1 and a = 1 for various grid sizes. 
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(c) N = 1000; a' = 0.3. 



(d) N = 1000; a 1 = 0.4. 
Figure 3. Concluded. 
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(d) Contour plot of figure 5(c). 
Figure 5. Continued. 
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(b) Contour plot at levels of 0.5, 1.0, and 1.5. 

Figure 7. Values of / distribution for (3 = 0.25. 
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(a) Classical motion. (b) Exact quantum motion. 



x 


(c) Monte Carlo simulation. 

Figure 8. Contour plots of fw/0 at levels of 0.5, 1.0, and 1.5 for 0 = 1, k = 1, Grid size = 0.3, and 
ol = 0.3 at t — 7 r for a pure state. The symbols x and p denote the position and momentum 
coordinates, respectively. 



(a) Classical motion. 


(b) Exact quantum motion. 



(c) Monte Carlo simulation. 

Figure 9. Contour plots of f u jf3 at levels of 0.5, 1.0, and 1.5 for /? = 1, k = 1, Grid size = 0.3, 
and a ' — 0.3 at t = 2n for a pure state. The symbols x and p denote the position and 
momentum coordinates, respectively. 




(a) Classical motion. 


(b) Exact quantum motion. 



(c) Monte Carlo simulation. 

Figure 10. Contour plots of fu>/0 at levels of 0.5, 1.0, and 1.5 for /3 = 1, k = 0.5, Grid size — 0.3, 
and a 9 = 0.3 at t — n for a weaker potential. The symbols x and p denote the position and 
momentum coordinates, respectively. 
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(a) Classical motion. 


(b) Exact quantum motion. 



(c) Monte Carlo simulation. 

Figure 12. Contour plots of fu/fl at levels of 0.5, 1.0, and 1.5 for /? = 0.5, 
Grid size = 0.3, and a' = 0.3 at t = 3 t\ for a mixed state. The symbols x and 
the position and momentum coordinates, respectively. 


k = 0.5, 
) denote 




(a) Classical motion. 


(b) Exact quantum motion. 


P 



(c) Monte Carlo simulation. 

Figure 13. Contour plots of f^jP at levels of 0.5, 1.0, and 1.5 for j3 — 0.5, k — 0.5, 
Grid size = 0.3, and o/ — 0.3 at t = 47t for a mixed state. The symbols x and p denote 
the position and momentum coordinates, respectively. 
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(c) Monte Carlo simulation. 

Figure 14. Contour plots of fuj/P at. levels of 0.5, 1.0, and 1.5 for (3 = 0.25, k — 0.5, 
Grid size = 0.3, and a ' = 0.3 at t = 3n for a mixed state. The symbols x and p denote 
the position and momentum coordinates, respectively. 


Approximate quantum motion 
Classical motion 



<x> 


Figure 15. First moments of Wigner distribution function (<x>, <p>) with k — 0.5 and j3 = 0.5 
over a period of time 0-47T. The cross (+) on the curves follows equal intervals of time. 
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